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ABSTRACT 


A  new  iterative  algorithm  for  the  maximum  entropy  power 
spectrum  estimation  is  presented  in  this  report.  The  algorithm 
which  is  applicable  to  two-dimensional  signals  as  well  as  one¬ 
dimensional  signals,  utilizes  the  computational  efficiency  of 
the  Fast  Fourier  Transform  (FFT)  algorithm  and  has  been 
empirically  observed  to  solve  the  maximum  entropy  power  spectrum 
estimation  problem.  Examples  are  shown  to  illustrate  the 
performance  of  the  new  algorithm. 
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I.  INTRODUCTION 

The  problem  of  power  spectrum  estimation  (PSE)  arises  in 
various  fields  such  as  speech  processing  [1] ,  seismic  signal 
processing  [2],  image  restoration  [3],  radar  [4],  sonar  [5], 
radio  astronomy,  etc.,  and  its  applications  range  from 
identifying  signal  source  parameters  and  transmission  channel 
characteristics  to  removing  noise  from  images  [3] .  Consequently, 
this  problem  has  received  considerable  attention  in  the 
literature  and  a  variety  of  techniques  for  power  spectrum 
estimation  have  been  developed. 

One  technique  which  has  been  studied  extensively  due  to  its 
high  resolution  characteristics  is  the  Maximum  Entropy  (ME) 
method.  For  one-dimensional  (1-D)  signals,  this  method  is 
equivalent  [7]  to  auto-regressive  signal  modelling,  and  thus 
it  leads  to  a  linear  problem  formulation  that  is  theoretically 
tractable  and  computationally  attractive  [7] .  Unlike  most 
other  PSE  techniques  such  as  the  conventional  methods  [8,9]  and 
the  Maximum  Likelihood  Method  [10],  however,  the  ME  method  does 
not  extend  from  1-D  signals  to  two-dimensional  (2-D)  signals  in 
a  straight-forward  manner  and  the  ME  method  for  2-D  signal  PSE 
remains  a  highly  non-linear  problem. 

To  solve  the  non-linear  ME  PSE  problem  for  2-D  signals, 
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various  attempts  [2,11,12,13]  have  been  made  in  the  literature. 
In  all  cases,  however,  the  algorithms  are  computationally 
unattractive,  and  there  is  no  guarantee  of  a  solution  or  only 
an  approximate  solution  can  be  obtained.  For  example,  Burg  [11] 
has  proposed  an  iterative  solution  which  requires  the  inversion 
of  a  matrix  in  each  iteration  where  the  dimension  of  the  matrix 
is  in  the  order  of  the  number  of  the  given  auto-correlation 
points.  No  experimental  results  using  this  technique  have  yet 
been  reported.  As  another  example,  Wernecke  and  D'Addario  [12] 
have  proposed  a  scheme  in  which  an  attempt  is  made  to  numeri¬ 
cally  maximize  the  entropy.  The  maximization  is  done  by 
continuously  adjusting  the  power  spectrum  (PS)  estimate  and 
evaluating  the  expressions  for  the  entropy  and  its  gradient. 

The  procedure  is  computationally  expensive  and  is  not  guaranteed 
to  have  a  solution.  As  a  third  example.  Woods  [2]  expresses 
the  ME  PS  estimate  as  a  power  series  in  the  frequency  domain 
and  attempts  to  approximate  the  ME  PS  estimate  by  truncating 
the  power  series  expansion.  Even  though  such  an  approach  has 
some  computational  advantages  relative  to  others,  the  method  is 
restricted  to  the  class  of  signals  for  which  the  power  series 
expansion  is  possible.  Furthermore,  examples  have  been  found 
in  which  the  algorithm  does  not  converge  to  the  desired  ME  PS 
estimate . 
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In  this  paper,  we  develop  a  new  iterative  algorithm  which 


is  computationally  simple  due  to  its  utilization  of  the  Fast 
Fourier  Transform  (FFT)  algorithm  and  which  leads  to  the  true 
ME  power  spectrum  estimate  for  2-D  signals  as  well  as  1-D 
signals.  In  Section  II,  we  review  briefly  previous  results  on 
the  ME  PSE  for  2-D  signals.  In  Section  III,  we  develop  a  new 
algorithm  for  the  ME  PSE.  In  Section  IV,  we  illustrate  and 
discuss  the  performance  of  this  algorithm  by  way  of  various 
examples . 

II.  PREVIOUS  RESULTS  ON  ME  POWER  SPECTRUM  ESTIMATION  FOR  2-D 
SIGNALS 


In  this  section,  we  review  briefly  important  previous 
results  relevant  to  this  paper  on  the  ME  PSE  for  2-D  signals. 
In  reviewing  these  results,  we  use  the  following  notations: 


x(n1,n2)  : 

Wn2|: 

Vnl'n2l! 

Px  (l «2)  1 

Px(ail'a)2>  : 
A (n1,n2)  : 


a  2-D  random  signal  whose  power  spectrum  we 
wish  to  estimate. 

auto-correlation  function  of  x(n^,n2) 
an  estimate  of  Rx(n^,n2> 
power  spectrum  of  x(n^,n2) 
an  estimate  of  px(“i'u2^ 

auto-correlation  function  whose  power  spectrum 
is  l/px  (<*>;,_,  1^2 ) 
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.  r  Alairii  1 1  Xh 


J 


A 


:  a  set  of  points  (n^,n2)  for  which  Rx(n^,n2) 


is  known 

P  :  discrete  time  Fourier  transform 

P  ^  :  inverse  discrete  time  Fourier  transform 

With  the  above  notations,  the  ME  PSE  problem  can  be  stated 
as  follows: 


Px(“l'“2) 


1 


(3) 


l  l  X  (n.  ,n~)  •e~I1“lnl-e-^ai2n2 
(n1#n2)6A 

and  Rx (n^  »n2)  *  F  1[Px(“i»w2)]  for  (n^n^SA  (4) 

The  above  problem  statement  for  the  ME  PSE  applies,  with 
appropriate  dimensionality  changes,  to  all  signals  regardless  of 
their  dimensionality.  The  solutions  to  the  problem,  however, 
strongly  depend  on  the  signal  dimensionality.  For  1-D  signals, 
the  mean  square  error  minimization  of  the  prediction  filter 
based  on  auto-regressive  signal  modelling  requires  solving  a 
set  of  linear  equations  for  the  filter  coefficients  and  the 
power  spectrum  obtained  from  the  estimated  filter  coefficients 
is  identical  to  the  ME  PS  estimate.  For  2-D  signals,  this  is 
no  longer  the  case.  Specifically,  even  though  minimizing  the 
mean  square  error  of  the  auto-regressive  filter  still  requires 
solving  a  set  of  linear  equations,  the  power  spectrum  obtained 
from  the  estimated  filter  coefficients  is  not  the  ME  PS  estimate. 
The  reason  for  this  can  be  seen  by  examining  the  form  of  the 
normal  equations  for  the  filter  coefficients  in  the  auto¬ 
regressive  signal  modelling.  The  derivation  of  the  general 
form  of  the  normal  equations  for  2-D  signals  is  analogous  to 
that  for  1-D  signals  and  is  given  by 
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where  represents  the  auto-regressive  filter  coefficients  to 
be  estimated,  the  set  B  consists  of  all  points  where  the  filter 
mask  has  non-zero  values,  and  the  power  spectrum  obtained  from 


From  equation  (5) ,  for  any  non-trivial  choice  of  B,  that  is  if 
B  does  not  consist  of  a  set  of  collinear  points,  the  size  of 
independent  values  of  R  (n^,^)  required  to  solve  the  above 
set  of  equations  is  greater  than  the  size  of  the  filter  mask. 

For  example,  consider  the  filter  mask  shown  in  Figure  1(a)  in 
which  the  solid  dots  represent  the  range  for  which  a^^  is  non¬ 
zero.  In  Figure  1(b)  is  shown  the  size  of  independent  values  of 
Rx  ( ni ,  n2  ^  rec3uire(3  to  solve  for  a^j  in  Figure  1(a)  by  equation 
(5).  Clearly,  the  number  of  correlation  points  needed  is 
greater  than  the  number  of  filter  coefficients.  Since  the 
estimated  power  spectrum  given  by  equation  (6)  is  completely 
determined  by  the  filter  coefficients  alone,  it  does  not 
possess  enough  degrees  of  freedom  to  satisfy  equation  (4)  which 
is  required  for  the  ME  PS  estimate.  Due  to  this  difficulty,  a 
closed  form  solution  for  the  2-D  ME  PSE  problem  has  not  yet 


(a)  First  quad 
b)  Independent  a 
mal  equations  fo 


ant  auto-regressi' 
to-correlation  po 
the  mask  of  (a ) . 


been  found. 

In  the  absence  of  a  closed  form  solution,  it  is  important 
to  know  the  conditions  for  the  existence  and  uniqueness  of  the 
solution.  In  this  regard,  Woods  [2]  has  obtained  the 
theoretical  result  that  if  the  given  R^(n1,n2)  for  (n^jn.JGA  is 
a  part  of  some  positive  definite  correlation  function  (meaning 
that  its  Fourier  transform  is  positive  for  all  (uj  ,u>  )  )  ,  a 
solution  to  the  ME  PSE  problem  exists  and  is  unique.  In  general, 
it  is  difficult  [14]  to  determine  if  the  given  segment  of  the 
correlation  function  is  a  part  of  some  positive  definite 
correlation  function,  even  though  this  is  generally  the  case 
in  most  practical  problems.  In  this  report,  we  assume  that 
the  given  segment  of  the  correlation  function  indeed  forms  a 
part  of  some  positive  definite  correlation  function  so  that 
the  solution  to  the  ME  PSE  problem  exists  and  is  unique. 

III.  A  NEW  ITERATIVE  ALGORITHM 

In  this  section,  we  develop  a  new  iterative  algorithm  for 
the  ME  PS  estimates  which  is  applicable  to  both  1-D  and  2-D 
signals.  This  algorithm  is  computationally  simple  since  it 
utilizes  the  computational  efficiency  of  the  FFT  algorithm. 
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Suppose  we  are  given  Rx^n^»n2^  f°r  (n^»n2^eA  such  that 


is  a  segment  of  some  positive  definite  correlation 


function.  To  find  the  unique  ME  PS  estimate  px (“^'“2^ '  we 
express  a  power  spectrum  Py(w^,<u2)  as  follows: 


py(“l'w2) 


=  F[Ry(n1,n2) ] 


00  00 


I  I  Rv(ni'n2) 

n^--°°  n2=-»  * 


-jto^n^  —  3w2^2 


and 


p  7 - r  —  F  [X  (n,  ,n»)  ] 

R  Ill)  .  rlil  t  I  J-  Z 


■  y  mu  1  r  uj  2 


<30  OO 


l  l  x(n  ,n2)* 

ni"-  n2=““ 


~^ini  ~3u2n2 
e  •  e 


(7) 


(8) 


From  equations  (7)  and  (8) ,  it  is  clear  that  R  (n.  ,  n~)  can 

y  i  z 

be  obtained  from  x  (n^n^  and  vice  versa  through  Fourier  trans- 

from  operations.  Now,  from  equations  (3)  and  (4)  Py (to ^ » u> 2 )  is 

the  unique  ME  PS  estimate  if  and  only  if  X(n^,n2)  =  0  for 

(n^n^gA  and  RyCn^n^  =  R^n^n^  for  (n^,n2)6A.  Thus,  we 

see  that  for  P  (<»>..  ,w_)  to  be  the  desired  ME  PS  estimate,  we 
y  1  z 

have  a  constraint  on  Ry(n^,n2)  and  a  constraint  on 
Recognizing  this,  it  is  straight-forward  to  develop  a  simple 
iterative  algorithm  to  find  the  unique  ME  PS  estimate. 
Specifically,  we  go  back  and  forth  between  Ry(n^,n2)  (the 
correlation  domain)  and  x (n, ,n„)  (the  coefficient  domain)  and 


at  each  time  impose  the  constraints  on  R  (n^,^)  and  X(n^,n2). 
Thus,  starting  with  some  initial  estimate  for  X(n^,n2)  we  obtain 
an  estimate  for  R^(n^,n2).  This  estimate  is  then  corrected  by 
the  given  Rx(n-^,n2)  over  the  region  A  and  is  used  to  generate 
a  new  x(n^,n2).  The  new  X(n^,n2)  is  then  truncated  to  the 
desired  limits  and  this  procedure  is  repeated.  The  above 
iterative  procedure  is  illustrated  in  Figure  2  and  forms  the 
basis  for  a  new  iterative  algorithm  for  the  ME  PSE. 

The  iterative  procedure  discussed  above  is  very  similar 
in  form  to  other  iterative  techniques  [15,16]  that  have  been 
successfully  used  in  image  processing.  Even  though  the'  condi¬ 
tions  under  which  the  algorithm  converges  are  not  yet  known,"5 
if  the  algorithm  converges  the  converging  solution  satisifies 
both  equations  (3)  and  (4)  and  consequently  is  tl}.e  desired  ME 
PS  estimate. 

The  algorithm  in  Figure  2  can  not,  in  general,  be  used  to 
obtain  the  ME  PS  estimate  without  some  modifications  due  to 
the  spectral  zero  crossing  problem.  Specif icilly,  the 
algorithm  in  Figure  2  requires  two  inversions  of  the  spectral 
estimates  in  each  iteration,  and  thus  the  algorithm  can  not  be 
continued  if  the  power  spectrum  estimate  has  a  zero  crossing  at 
any  stage  in  the  iterative  procedure.  Unfortunately,  zero 
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1 02455-N 


INITIAL  ESTIMATE  OF  XfrynJ 


Ry(ni  ,n2>  =  F 


{f  [X(n1,n2))} 


I 


CORRECT  R  (n1#n2)  WITH  R^ryr^)  FOR  (n^) 

X(nrn2)  =  F  1{pTRy(vn2)]} 

I 

X^^tij)  =  0  FOR  (n^n^/A 

Px^uru2)  =  F  !Ry(ni'n2)] 


Fig.  2.  A  new  approach  to  2-D  maximum  entropy 
(ME)  power  spectrum  estimation  (PSE) . 
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crossings  can  occur  in  two  different  ways  in  each  iteration. 
One  is  the  correction  of  the  correlation  function  and  the 
other  is  the  truncation  of  the  coefficients.  To  see  this,  let 
Xm(n1,n2)  and  R^n^n^  represent  Mn^n^  and  R  (n^rij)  after 
mth  iteration  and  suppose  that  the  following  conditions  hold; 


F  [xm(n1,n2)  ]  >  0  for  all  (w^u^). 

(9) 

F  [R™  (n^  ,n2 )  ]  >  0  for  all  (o)^,<»>2), 

(10) 

and  xm(n,,n,)  =f"1[ - - — - - ]-w(n,  ,n,) 

1  2  F[R^(ni,n2)] 

(ID 

where  w(n1#n2)  represents  a  rectangular  type  window  such  that 


w (n^ ,n2)  =  1  for  (n^n^SA 
0  otherwise 


(12) 


Similarly,  let  xm+1  (n^n^  and  R^+1(n1(n2)  represent  X  (nx ,  n2) 

and  R  (n, ,nj  after  m+lth  iteration.  In  the  iterative  algorithm 
y  1  2 

of  Figure  2,  xm+1(n1,n2)  and  R^+1(n1,n2)  are  obtained  from 


Xm(nx,n2)  by 


R' (n1,n2)  -  F  ^ [ 


F [Xm(n1,n2) ] 


(13) 


Rm+1(n1,n2)  =  Rx(n1'n2)  for  (n1»n2)eA 


R' (n^ ,n2)  otherwise 


(15) 


A  1  (n,  ,n_)  =  F-1! - jrrr^ - )/ 

1  2  F[R^+i(n1,n2)] 

and  Xm+1(n1<n2)  =  X'(n1,n2)  for  (n^n^GA 

0  otherwise 

=  X ' (n1,n2) -w (n1,n2)  (16) 

From  equations  (13)- (16),  it  is  clear  that  R'(n^,n2)  is 
positive  definite  since  Am(n^,n2)  is  assumed  to  be  positive 

_  i  1 

definite  but  R  (n, ,n«)  may  not  be  positive  definite  due  to 

y  -*•  *■ 

the  rectangular  windowing  w(n^,n2)  in  equation  (14).  Further¬ 
more,  even  if  Rm+1(n..,n  )  were  positive  definite  so  that 

y  i  * 

m  i  "I 

X'(n1,n2)  is  positive  definite,  A  (n^n^  may  not  be 
positive  definite  due  to  wCn^n^)  in  equation  (16). 

To  ensure  the  resulting  ft™*1  (n^  ,n2)  and  \m+^(n^,n2)  are 
positive  definite  so  that  the  iterations  can  be  continued,  we 
make  a  modification  to  equation  (14)  .  Specifically,  instead  cf 
forming  R™+^ (n^,n2)  by  replacing  R'(n^,n2)  with  all  its  known 
values,  namely  Rx(n]/n2^  for  (n^,n2)6A  in  equation  (14),  suppose 
we  form  Ry+^  (n^,n2>  by  linearly  interpolating  between  the  values 
of  R'(n^,n2)  and  the  known  values  of  Rx(n^,n2)  f°r  (n^,n2)€A. 
Then,  in  the  modified  iterative  algorithm,  A™*'*' (n^  ,n2)  and 
Ry+1(nlfn2)  are  obtained  from  Xm(n^,n2)  by 


R’(n  ,n  )  =  F_1[ - ^ - 

1  Z  F[Am(n1,n2)) 
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and 


Rym+1(nl'n2) 


X ' (ni'n2* 


,  m+1 ,  . 

X  (n.,n  ) 


=  a  *R’  (n1#n2)  +  (1-a  )  *R  ,n2)  for  (n^n^eA 
R1 (n^ , n2 )  otherwise 
=  R’ (nx ,n2) + (1-a ) ' (R^ (n^ , n2 ) -R 1 (n1 ,n2) ) 


•w (n1 ,n2) 

(18) 

-1  1 

F(Rym+1(n1,n2)] 

(19) 

=  X ' (n1 , n2 ) -w(n1#n2) 

(20) 

Comparing  equations  (14)  and  (18) ,  equation  (18)  reduces  to 
equation  (14)  when  a  =  0.  With  any  other  choice  of  a,  equation 
(18)  represents  a  non-ideal  correction  of  R' (n^,n2)  with  the 
known  values  Rx(ni»n2^  f°r  )6A,  with  a  larger  deviation 

of  a  from  zero  corresponding  to  a  more  non-ideal  correction. 
However,  with  proper  choice  of  a,  the  resulting  R  m+1 (n. , n_ ) 

y  x  z 

and  Xm+  (ni,n2)  can  be  guaranteed  to  be  positive  definite. 

This  can  be  seen  by  noting  that  Xm(n^,n2)  and  therefore 
R'(n^,n2)  are  assumed  to  be  positive  definite  and  by  considering 
a  sufficiently  close  to  1  so  that  R^.m+'*' (n^^  ,n2)  and  Xm+1(n^,n2) 
can  be  made  arbitrarily  close  to  R'(nlfn2)  and  Xm(nlfn2). 
Therefore  by  properly  choosing  a  in  the  range  0  <_  a  <  1,  the 
spectral  zero  crossing  problem  can  be  avoided  and  the  iterations 
can  be  continued. 
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From  equations  (13)— (16) ,  it  is  clear  that  if  Xm(n^,n2>  and 

R-(ni,n2)  satisfy  equations  (9)— (11) ,  then  Xm+^(n^,n2)  and 
+  1 

(n^,n2)  obtained  by  the  modified  iterative  algorithm  also 
satisfy  equations  (9) -(11).  With  proper  choice  of  a,  then, 
if  X°(n^,n2)  and  R^(n^,n2),  the  initial  estimates  of  X(n^,n2) 
and  Ry(n^,n2),  satisfy  equations  (9)~(11),  the  iterations 
specified  by  equations  (9) -(12)  and  (17) -(20)  form  an 
iterative  algorithm.  This  algorithm  is  shown  in  Figure  3. 

In  implementing  the  algorithm  in  Figure  3,  there  are 
several  important  issues  that  need  to  be  discussed.  One  of 
them  is  the  determination  of  the  length  of  the  Discrete  Fourier 
Transform  (DFT)  and  Inverse  Discrete  Fourier  Transform  ( IDFT)  to 
be  used  for  the  Fourier  transform  and  inverse  Fourier  transform 
operations.  In  general,  a  large  DFT  length  should  be  used  in 
the  implementation  to  avoid  any  aliasing  problem.  Specifically, 
the  ME  method  of  PSE  is  essentially  an  attempt  to  extrapolate 
the  correlation  function  beyond  the  limits  of  the  known  segment. 
Since  the  DFT  is  used  in  the  implementation  instead  of  the  true 
Fourier  transform,  the  length  of  the  DFT  should  be  chosen  such 
that  the  extended  correlation  function  corresponding  to  the 
ME  PS  estimate  is  essentially  zero  beyond  the  DFT  limits. 

Choice  of  the  DFT  length  and  its  effect  on  the  system  perform¬ 
ance  will  be  further  discussed  in  Section  IV. 
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Another  issue  to  be  considered  in  the  implementation  is 
determination  of  and  Ry(n^,n2),  the  initial  estimates 

of  X(nlfn2)  and  R^n^n^.  Even  though  various  different 
choices  are  possible,  one  choice  which  is  particularly  simple 
and  satisfies  equations  (9) -(11)  is 

Ry^l'1^)  =  Rx(0'0)  *6  (ni'n2)  (21) 

and  X°(n1,n2)  =  R  (o  ~  o')~  ’ 6  (ni  >n2  *  (22) 

where  it  is  assumed  that  the  region  A  includes  the  point 
^nl,n2^  =  an^  Rx(0,0)  >  0.  Furthermore,  the  choice  of 

X°(n1,n2)  and  Ry(n1,n2)  given  by  equations  (21)  and  (22)  is 
reasonable  in  that  R™^,^)  =  *w(n^,n2)  for  m  =  1  and 

a  =  0  in  the  iterative  algorithm  of  Figure  3. 

A  third  issue  to  be  discussed  in  the  implementation  is  how 
the  value  of  a  is  specifically  determined  in  each  iteration. 

The  choice  of  a  is  dictated  by  two  considerations.  One  is  the 
requirement  that  the  resulting  R^m+1 (n^,n2)  and  Xm+^(n^,n2)  are 
positive  definite.  The  second  is  our  desire  to  choose  a  as 
close  to  zero  as  possible  so  that  more  correction  with  the 
known  correlation  points  is  made  in  each  iteration.  Thus,  the 
ideal  choice  of  a  is  the  smallest  a  in  the  range  0  <_  a  <  1 
for  which  the  resulting  Rym+1  (n^.nj)  and  Xtn+1(n1,n2)  are 
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positive  definite.  Finding  the  optimum  value  of  a,  however,  may 
require  many  iterations  in  which  we  begin  with  a  =  0  and 
successively  increase  a  until  the  resulting  Ry.m+^(n^,n2)  and 
Xm+^ (n^ ,n2)  are  positive  definite.  Since  each  iteration  requires 
one  DFT  and  one  IDFT ,  obtaining  a  alone  can  be  a  significant 
computational  burden.  Further,  it  has  been  empirically 
determined  that  choosing  the  smallest  a  in  each  iteration  can 
lead  to  a  limit  cycle  behavior  where  the  algorithm  does  not 
converge  to  the  ME  PS  estimate. 

An  alternative  approach  to  the  choice  of  a ,  which  avoids 
the  above  two  problems  and  is  used  in  this  paper,  is  to  begin 
with  0^  =  0  and  change  a  in  the  following  manner; 

min  [F  [R’  (n.  ,n  )  ]  ]  (23) 

( ,  w2 ) 

“m+l^^m'1'1''  I'  min  F  t  (R  ,  n  )  -R1  (n.,n2)  )  *w(n,,n  )  ]  I  ] 

1 ,w2)  1 

if  min  [F [ (R^ ^ ,n  ) -R' ^ , n2) ) • w (n  , n2) ) ]  <  0 
( u)  ^  ,  co  2 ) 

and  a  ,,*-a  otherwise, 
m+1  m 

where  represents  the  value  of  a  in  the  ith  iteration, 

max  [  ,  ]  represents  the  maximum  of  the  two  arguments,  min  [  ] 

(oj^,w2) 

represents  the  minimum  of  the  argument  expression  over  (ui^,oj2), 
and  "k"  is  a  scalar  constant  with  0<k£l.  Typically,  "k"  is 
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chosen  to  be  0.01,  but  can  be  increased  to  1  for  low  S/N  ratio 
cases.  When  a  is  chosen  according  to  equation  (23),  it  is 
straight-forward  to  show  from  equation  (14)  that  the  resulting 
R^ra+^ (n^ , n^ )  is  always  positive  definite.  In  addition,  it  has 
been  empirically  observed  that  the  resulting  xm+^(n^,n2)  is 
almost  always  positive  definite.  In  those  rare  cases  when 
Xm+^ (n^ , n^ )  is  not  positive  definite,  a  is  further  increased 
and  the  iterations  continue  with  the  increased  a . 


Finally,  another  important  issue  to  be  considered  in 
implementing  the  algorithm  in  Figure  3  is  the  decision  on  when 
the  algorithm  converged  so  that  the  iteration  can  be  stopped. 
One  reasonable  approach  is  to  consider  that  the  algorithm 
has  converged  when  the  following  condition  is  satisfied. 


I  I  (R*  (n  ,n  ) -R  (n.  ,n_)  )  2 

(nlfn2)eA  1  ^  X  l 


l  l  4 

(ni,n2)6A  x 


(n^,n2) 


< 


e 


(24) 


clearly,  if  e  =  0  and  the  algorithm  has  converged,  then  the 
converging  solution  corresponds  to  the  desired  ME  PS  estimate. 
However,  due  to  a  finite  DFT  length  and  finite  precision 
arthematic  used,  it  is  not  possible  to  reduce  the  error 


exactly  to  zero.  A  convergence  decision  is  made  when  the  error 
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e  is  very  small,  typically  10 

In  Figure  4  is  shown  a  more  detailed  flow  chart  of  an 
algorithm  which  incorporates  the  important  implementation 
issues  discussed  above.  It  is  not  theoretically  known  under 
what  conditions  the  algorithm  in  Figure  4  converges.  As  will 
be  discussed  in  the  next  section,  however,  we  have  empirically 
observed  that  the  algorithm  always  converges  to  the  ME  PS 
estimate  with  a  sufficiently  large  choice  of  the  DFT  length  and 
a  sufficiently  small  choice  of  e. 

IV.  EXAMPLES  AND  DISCUSSIONS 


In  this  section,  we  illustrate  and  discuss  various 
examples  in  which  the  algorithm  in  Figure  4  has  been  applied 
to  obtain  the  ME  PS  estimate.  In  all  cases,  it  was  assumed 
that  the  correlation  function  originated  from  sinusoids  buried 
in  white  noise  so  that  the  correlation  function  given  has  the 
form  of 


M 


2  **  2 

R  (n)  =  a  • 6 (n) +  V  a • • cos (w . n)  for  1-D  signals  (25) 

x  i=l  1 


M 


and  Rx  (n^ ,  n2 )  =  a  •  S  (n^n^  +  £  aScostio.^n^u.^n^ 


for  2-D  signals 


(26) 


where 


2 

o 


represents  the  white  noise  power,  M  represents  the 
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Fig . 
for  2 


DFT  LENGTH  ’N',  SCALE  FACTOR  ’K  ’ 

a  --  0.0,  «  10'4 

1 

INITIAL  ESTIMATE:  ▼ 


R°(n1,n2)  =■  Rx(0,0)  -  6(n,,n2) 


*71  ,  v 

R*(n  n9)  =  i  DFT  « - - -  V 


Loft  i> 


DEFINE:  • 

=  min  { DFT  j 


r  2  i 

X2  =  min  (DFT  ( (R^tn^  ,n2>  -  R'fn^n^]-  w(n^,n2)|) 


S'"!1 

YES 


* 


NO 

nrn'a  \l  -iTT*' 


"I 


R  (n^nj)  =  R'(n],n2)  *  (I  -  °mt]  )  iR^",  ,n2)  -  R'fn,,^))  *  w(n},n2> 


V  i 


X™  ^ (n1  /n2>  =  '  w(nj/n2^ 

DFT[Xm*,(n)#n2)t  >  0  ?  — ° - 


J 


YESj 

^  1  |R,<nl'n2)  "  R'<Vn2!|2 
nl  n2 


Px(Ul,a2)  =  OfT  fR  Vy nj)l 

1.  A  detailed  flow-chart  of  the  iterative  algorithm 
-D  ME  PSE  implemented  in  this  paper. 
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number  of  sinusoids,  represents  the  power  of  the  ith 
sinusoid,  and  for  1-D  signals  and  and  for  2-D 

signals  represent  the  frequencies  of  the  ith  sinusoid. 

We  first  consider  the  case  of  1-D  signals.  For  1-D 
signals,  the  iterative  algorithm  in  Figure  4  is  not  very  useful 
in  obtaining  the  ME  PS  estimate  due  to  the  existence  of  a 
simple  closed  form  solution.  However,  the  1-D  signal  case  is 
ideal  in  illustrating  that  the  solution  obtained  from  the 
iterative  algorithm  in  Figure  4  indeed  leads  to  the  ME  PS 
estimate.  In  Figure  5  are  shown  the  results  obtained  using  the 
parameters  in  Table  1.  In  Figure  5(a)  are  shown  the  results 
obtained  from  the  iterative  algorithm  as  a  function  of  the 
number  of  iterations.  The  converging  solution  with  the 
choice  of  e  in  Table  1  was  obtained  after  200  iterations.  In 
Figure  5 (b)  is  shown  the  ME  PS  estimate  obtained  from  the 
closed  form  solution  of  equation  (6).  Figure  5(c)  shows  the 
PS  estimates  obtained  from  the  iterative  algorithm  (solid  line) 
and  the  closed  form  solution  (dotted  line) .  It  is  clear  from 
Figure  5(c)  that  the  iterative  algorithm  leads  to  the  ME  PS 
estimate . 

As  another  1-D  example,  Figure  6  is  similar  to  Figure 
5  except  that  the  parameters  in  T’able  2  were  used  to  generate 
the  results.  In  addition  to  the  above  two  examples,  a  variety 
of  other  examples  have  been  considered.  In  all  cases 


22 


TABLE  1 


Parameters  of  a  one-dimensional  example  for  Figure  5  . 


M 

2 

a 

2 

ai 

w-^/2-n 

e 

NDFT 

NITR 

TIME 

1 

1.0 

1.0 

0.1 

0.3456 

10-4 

128 

200 

3  secs. 

A  :  size  of  the  known  auto-correlation  array, 
symmetric  about  the  origin 

M  :  number  of  sinusoids 
2 

o  :  noise  power 
2 

i  :  power  of  ith  sinusoid 

<iK  :  frequency  of  ith  sinusoid 

e  :  error  used  for  the  convergence  test 

NDFT  :  size  of  discrete  Fourier  transform  used 

NITR  :  number  of  iterations  required  to  reach  e 

Time  :  the  CPU  time  required  using  IBM-370  at  M.I.T. 
Lincoln  Laboratory 
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0.0  0.1  0.2  0.3  0.4 

ID/  2 TT 

(a)  The  PS  estimate  as  a  function  of  the  number  of 
iterations  (NITR) .  NITR  =  0,  25,  50,  75,  100,  125, 
150. 


Fig.  5.  1-D  PSE  for  the  data  in  Table  1. 
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0.1  0.2  0.3  0.4  ,,  .  0.5 

u)/2n  -*■ 

(b)  Results  of  ME  PSE  by  the  closed  form  solution 
of  equations  (5)  and  (6)  . 


Fig.  5.  Continued. 
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(c)  Result  of  the  iterative  method  (solid  line) 
Result  of  the  closed  form  solution  (dotted  line) 
The  two  results  are  indistinguishable. 


Fig.  5.  Continued. 


that  we  have  considered  so  far,  the  iterative  algorithm  in 
Figure  4  lead  to  the  PS  estimates  which  are  visually 
indistinguishable  from  the  closed  form  solutions. 


TABLE  2 

Parameters  of  a  one-dimensional  example  for  Figure  6. 
Notations  used  in  the  table  are  explained  in  Table  1. 


D 

M 

_ 

2 

"1 

2 

ai 

. 

_ 

NDFT 

NITR 

Time 

5 

D 

1 

10 

0.1 

mm 

64 

69 

less  than  1  sec. 

We  now  consider  the  case  of  2-D  signals.  For  2-D  signals, 
a  closed  form  solution  for  ME  P3E  has  not  yet  been  found  and 
consequently  the  iterative  algorithm  developed  in  this  paper 
has  practical  significance.  In  Figure  7  are  shown  the  results 
obtained  using  the  parameters  in  Table  3.  In  Figure  7(a)  is 
shown  the  periodogram  [9]  obtained  by  Fourier  transforming 
Rx(Vn2)-w(Vn2>.  In  Figure  7(b)  is  shown  the  PS  estimate 
based  on  auto-regressive  signal  modelling  of  equation  (6)  using 
a  backward  L  shaped  filter  mask  [16]  .  This  particular  mask  was 
chosen  due  to  its  high  performance  [16]  relative  to  other  filter 
mask  shapes  such  as  symmetric  and  first  quadrant  filters.  In 
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0.0 

(a) 

NITR 

Fig 


J - 1 - 1  ■  - 

0.1  0.2  0.3  0.4  ..  0.5 

u)/2u  ■+■ 


The  PS  estimate  as  a  function  of  NITR. 
=  12,  24,  36,  48,  60. 


6.  1-D  PSE  for  the  data  in  Table  2. 
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(b)  Result  of  ME  PSE  by  the  closed  form  solution 
of  equations  (5)  and  (6)  . 


Fig.  6.  Continued. 


(c)  Result  of  the  iterative  method  (solid  line) 
Result  of  the  closed  form  solution  (dotted  line) 
The  two  results  are  indistinguishable. 


Fig.  6.  Continued. 


Figure  7(c)  is  shown  the  PS  estimate  obtained  from  the 
iterative  algorithm  developed  in  this  paper.  To  ensure  that 
the  PS  estimate  shown  in  Figure  7 (c)  corresponds  to  the  ME  PS 
estimate,  an  additional  estimate  was  obtained  using  a  larger 

__  g 

DFT  size  of  128  *  128  points  and  a  much  smaller  e  of  10  and 
was  compared  to  the  result  in  Figure  7(c).  The  two  estimates 
were  visually  indistinguishable.  It  is  clear  from  Figure  7 
that  the  two  sinusoids  are  resolved  only  in  Figure  7(c),  implying 
that  the  ME  method  for  PSE  has  the  high  resolution  characteris¬ 
tics  for  2-D  signals  as  well  as  1-D  signals. 


TABLE  3 

Parameters  of  a  two-dimensional  example  for  Figure  7. 
Notations  used  in  the  table  are  explained  in  Table  1. 
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In  using  the  iterative  algorithm  developed  in  this  paper, 

care  should  be  exercised  in  properly  choosing  the  DFT  length. 

As  has  been  discussed  in  Section  III,  a  smaller  aliasing  error 

requires  a  longer  DFT.  However,  once  the  DFT  length  is  chosen 

so  that  the  aliasing  error  is  sufficiently  small  relative  to  c ,  a 

further  increase  in  the  DFT  length  makes  little  improvement  in 

the  spectral  estimates,  but  only  increases  the  computational 

requirements.  To  give  a  rough  estimate  of  the  DFT  length  needed 

in  practice,  a  number  of  examples  have  been  considered  and  the 

-4 

ranges  of  DFT  lengths  required  for  e  =  10  have  been  computed 
as  a  function  of  the  S/N  ratio  and  the  size  of  A.  The  S/N 
ratio  in  dB  is  defined  as 

M  2 

-l  Kl 

S/N  ratio  in  dB  A  10«log(- - j - )  (27) 

and  the  computed  results  are  shown  in  Table  4. 


TABLE  4 

Typical  ranges  of  minimum  DFT  size  needed  in  practice 
to  achieve  convergence  of  the  algorithm  in  Figure  4 
with  e  =  10 


Size  of  A 


S/N  ratio 

1-D 

Signals 

2-D 

Signals 

5 

9 

3x3 

5x5 

-10  db 

8-16 

8-32 

8x8 

8x8  -  16x16 

0  db 

16-64 

32-128 

8x8  -  16x16 

16x16  -  64x64 

+10db 

64-512 

64-1024 

32x32-256x256 

64x64-512x512 

mammm 


From  Table  4,  it  is  clear  that  cases  in  which  the  S/N  ratio 

is  higher  or  the  size  of  A  is  larger  require  longer  DFTs .  This 

is  due  to  the  fact  that  the  correlation  functions  given  in 

equations  (25)  and  (26)  represent  more  sinusoidal  behavior  in 

region  A  for  a  higher  S/N  ratio  or  larger  size  of  A.  Therefore, 

the  extended  correlation  functions  by  the  ME  method  for  such 

cases  tend  to  be  longer  and  thus  require  longer  DFTs.  To 

further  illustrate  this  point  with  specific  examples,  we  have 

2 

considered  a  case  in  which  a  =0.5  with  all  other  parameters  to 
be  the  same  as  in  Table  3  and  another  case  in  which  the  size  of 
A  is  7x7  with  all  other  parameters  to  be  the  same  as  in  Table  3. 
To  obtain  the  desired  ME  PS  estimate,  both  cases  required  the 
DFT  size  of  128x128. 

If  there  is  significant  aliasing  error  relative  to  the 
e  used  for  the  convergence  test,  then  the  algorithm  does  not 
converge.  For  example,  if  the  DFT  size  of  32x32  is  used  for 
the  example  of  Figure  7,  the  algorithm  does  not  converge.  In 
such  a  case,  the  DFT  length  has  to  be  increased  to  obtain  the 
ME  PS  estimate.  If  the  DFT  length  can  not  be  increased  due 
to  computational  constraints,  e  has  to  be  increased  to  toler¬ 
ate  a  larger  aliasing  error.  In  such  a  case,  the  resulting 
PS  estimate  will  only  be  an  approximation  to  the  ME  PS 
estimate . 
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In  general,  cases  in  which  the  S/N  ratio  is  lower  or  the 

size  of  A  is  smaller  require  less  computations  to  reach  the 

desired  ME  PS  estimate  for  two  reasons.  One  reason  is  that  such 

cases  require  shorter  DFT  lengths,  as  are  clear  from  Table  4. 

2 

Another  reason  is  that  a  larger  noise  power  (o  )  contributes 
a  larger  positive  spectral  component  in  the  correlation  correct¬ 
ion  step  and  the  Fourier  transform  of  a  smaller  rectangular 
window  has  a  smaller  amplitude  in  its  negative  lobe.  The  effdct 
of  this  is  generally  a  smaller  value  of  a  in  each  iteration  and 
therefore  more  ideal  correction  of  the  correlation  function  in 
each  iteration.  Consequently,  for  such  cases,  a  smaller 
number  of  iterations  is  required  to  reach  the  ME  PS  estimate. 

In  fact,  when  the  S/N  ratio  is  sufficiently  low  and  the  size  of 
A  is  sufficiently  small,  the  value  of  a  can  be  chosen  to  be  0 
in  all  iterations  and  the  computation  time  can  be  significantly 

reduced.  As  an  example,  we  have  considered  a  case  in  which 
2 

o  =6.0  and  size  of  A  =  3x3  with  all  other  parameters  to  be 
the  same  as  in  Table  3.  The  computation  time  required  in 
generating  the  ME  PS  estimate  for  this  case  was  4.5  seconds  of 
CPU  time  using  IBM  370  at  M.I.T.  Lincoln  Laboratory.  This  is 
significantly  shorter  than  75  seconds  required  to  generate  the 
result  in  Figure  7(c). 
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Two  additional  examples  of  2-D  ME  PS  estimates  are 
shown  in  Figures  8  and  9.  The  results  in  these  two  figures 
were  obtained  using  the  parameters  in  Tables  5  and  6  respec¬ 
tively  . 

In  addition  to  the  above  examples,  we  have  considered  a 
variety  of  others.  In  all  cases  that  we  have  considered  so 
far,  we  have  empirically  observed  that  the  algorithm  always 
converges  to  the  ME  PS  estimate  for  a  sufficiently  large 
choice  of  DFT  length  and  a  sufficiently  small  choice  of  e. 

In  addition,  we  have  also  observed  that  the  ME  PS  estimates 
for  2-D  signals  have  high  resolution  characteristics  similar 
to  the  1-D  case.  A  more  quantitative  study  on  the  high 
resolution  characteristics  of  the  ME  PSE  method  for  2-D 
signals  is  currently  under  investigation. 

In  this  paper,  we  have  developed  a  specific  numerical 
algorithm  to  estimate  the  power  spectrum  of  a  signal  by  the 
ME  method.  Even  though  this  algorithm  led  to  successful 
results,  there  is  considerable  room  for  further  imporvements 
of  the  algorithm.  For  example,  to  avoid  the  spectral  zero 
crossing  problem  in  both  the  correlation  domain  and  the 
coefficient  domain,  we  have  considered  a  smaller  correction 
only  in  the  correlation  domain.  An  alternative  approach 


38 


TABLE  5 


Parameters  of  a  two-dimensional  example  for  Figure  8. 
Notations  used  in  the  table  are  explained  in  Table  1. 
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TABLE  6 

Parameters  of  a  two-dimensional  example  for  Figure  9. 
Notations  used  in  the  table  are  explained  in  Table  1. 
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would  be  to  perform  a  smaller  correlation  correction  and  a 
smaller  coefficient  correction  to  avoid  the  spectral  zero 
crossing  problems  in  the  correlation  and  coefficient  domains 
respectively.  Such  an  approach  may  allow  more  optimum  choice 
of  the  parameters  such  as  a  without  significant  computational 
burden.  This  and  other  ways  to  improve  the  performance  of 
the  algorithm  are  currently  under  investigation. 

Finally,  in  this  report,  we  have  considered  the  ME  PSE  of 
only  1-D  and  2-D  signals.  However,  the  iterative  algorithm 
developed  is  based  on  the  notion  that  the  ME  PS  estimate  should 
be  consistent  with  the  given  correlation  points  in  the  region 
A  and  the  corresponding  coefficients  should  be  zero  outside  the 
region  A,  which  car.  be  applied  to  signals  of  all  dimensions. 
Consequently,  the  iterative  algorithm  developed  in  this  paper 
may  also  be  useful  for  the  ME  PSE  of  signals  whose  dimensions 
are  higher  than  two. 
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